Dewetting-controlled binding of ligands to hydrophobic pockets 
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Abstract 

We report on a combined atomistic molecular dynamics simulation and implicit solvent analysis of 
a generic hydrophobic pocket-ligand (host-guest) system. The approaching ligand induces complex 
wetting/dewetting transitions in the weakly solvated pocket. The transitions lead to bimodal 
solvent fluctuations which govern magnitude and range of the pocket-ligand attraction. A recently 
developed implicit water model, based on the minimization of a geometric functional, captures 
the sensitive aqueous interface response to the concave-convex pocket-ligand configuration semi- 
quantitatively. 
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The water-mediated interaction between a ligand and a hydrophobic binding pocket 
plays a key role in biomolecular assembly processes, such as protein-ligand recogni- 
tion , the binding of the human immunodeficiency virus (HIV) [7] or the 
dengue virus |8j to human cells, the inhibition of influenza virus infectivity [9|, or in synthetic 
host-guest systems [10] . Experiments and explicit- water molecular dynamics (MD) simula- 
tions suggest that the concave nature of the host geometry imposes a strong hydrophobic 



constraint and can lead to very weakly hydrated pockets 
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to nanoscale capillary evaporation triggered by the approaching ligand 
called 'dewetting' transition has been also observed in other (protein) geometries 
It has been speculated that dewetting may lead to a fast host-guest recognition accelerating 
the hydrophobic collapse and binding rates of the ligand into its pocket . A deeper 

physical understanding of these sensitive hydration effects in hydrophobic recognition is still 
elusive. 

On the coarse-grained modeling side, the thermodynamics of molecular recognition is 
typically approached by surface area (SA) models [jjj]. A major flaw of these implicit sol- 
vent models is that the aqueous interface around the macromolecules is predefined (typically 
by rolling a probe sphere over the van der Waals surface) and is therefore a rigid object that 
cannot adjust to local energetic potentials and changes in spatial molecular arrangements. 
In particular, the dewetting 
trostatics, and geometry 11 . 
models. Their qualitative deficiency to describe the hydrophobic pocket- igand interaction 
in proteins [16], pocket models [13], or dewetting in protein folding [IT]] is therefore not 
surprising. 

In this letter, we combine explicit-water MD simulations and the variational implicit sol- 
vent model (VISM) [15] applied to a generic pocket-ligand model. The simulations show that 
the approaching ligand first slightly stabilizes the wet state in the weakly hydrated pocket, 
whereas, upon further approach, bimodal fluctuations in the water occupancy of the pocket 
are induced, followed by a complete pocket dewetting. The onset of fluctuations defines the 
critical range of pocket-ligand attraction. The VISM calculation, based on the minimization 
of a geometric functional, reproduces the bimodal hydration and explains it by the existence 
of distinct metastable states which correspond to topologically different water interfaces. As 
opposed to previous (SA type of) implicit models, VISM captures the range of the pocket- 



transition, which is highly sensitive to local dispersion, elec- 
, Ji]], can, per definitionem, not be captured by SA type of 
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ligand interaction semi-quantitatively. Strikingly, the observed nanoscale phenomena can 



be thus explained by geometric capillary effects, well-known on macroscales [18]. Explicit 
inclusion of dispersion interactions and curvature corrections, however, seem to be essential 
for an accurate description on nanoscales. 

Our generic pocket-ligand model consists of a hemispherical pocket embedded in a rect- 
angular wall composed of neutral Lennard- Jones (LJ) spheres interacting with U^(r) = 
4e[(o~/r) 12 — (o~/r) 6 ]. The atoms are aligned in a hexagonal closed packed arrangement with 
a lattice constant of 1.25 A. The LJ parameters are chosen to model a paraffin-like material 



and are e p = 0.03933 kJ/mol and a p = 4.1 A [13j, |27|]. We consider two different pocket 
radii: R = 5 and 8 A, which we refer to in the following as 'R5' and 'R8' systems. The 
ligand is taken as a methane (Me) represented by a neutral LJ sphere with parameters 
e = 0.294 kJ/mol and a = 3.730 A. It is placed at a fixed distance d from the flat part 
of the wall surface (z = 0), along the pocket symmetry axis in ^-direction, see Fig. 1 (a) 
for an illustration. The explicit-water MD simulations are carried out with the CHARMM 
package [3] employing the TIP4P water model in the NVT ensemble, two anti-symmetric 
walls with thickness of 7.5 A in a surface-to-surface distance of 30 A in a rectangular box 
with lengths L x = L y ~ 34 A and L z = 100 A, and 3D particle mesh Ewald summation. 
In equilibrating simulations, the volume V of the system was varied until the density in the 
center of the slab matched the bulk density of TIP4P water at a pressure of P = 1 bar 
and temperature T = 298 K. More technical details of the simulation setup can be found 
elsewhere [3,0]. A MD simulation snapshot is shown in Fig. 1 (b). 



The VISM was introduced in detail previously 15] and applied to the solvation of nonpolar 
solutes 2(]] . Briefly, let us define a subregion V void of solvent in total space W, for which we 
assign a volume exclusion function v(f) = for r*e V and v(f) = 1 else. The volume V and 
interface area S of V can then be expressed as functionals of v(r) via V[v] = f w d 3 r [1 — v(r)] 
and S[v] = J w d 3 r \Vv(r)\ = J dw dS, and the solvent density is p(f) = p v(r), where po is 
the bulk value. The solvation free energy G is defined as a functional of the geometry v(f) 



of the form 15] 



G[v] = PV[v]+ [ dS llv [l-25H(f)] 

JdW 



+ po d 3 r v{r)U{r), (1) 
Jw 

where ji v is the liquid-vapor interface tension, 5 the coefficient for the curvature correction 



of 7^ in mean curvature H(r), and U(r) = Ylf s ^^{r — fi) sums over the LJ interactions of 
all N s solute atoms at (ligand+wall atoms) with the water. The 5-term in (1) has been 
used in scaled-particle-theory 2l| for convex solutes only, generalized capillary theory 231 ] . 



and in the morphometric approach applied to the solvation of model proteins 
minimization SG[v]/5v = leads to the partial differential equation (PDE) 15] 



22j. The 



P - 2 7lv \H{r) + 8K{f)\ - Po U{r) = 



(2) 



which is a generalized Laplace equation of classical capillarity 18|, |23| extrapolated to mi- 
croscales by dispersion and the local Gaussian curvature K(r). The PDE (2) is solved using 
the level-set method which relaxes the functional (1) by evolving a 2D interface in 3D space 
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24j. The 
59 mJ/m 2 for 



and robustly describes topological changes, such as volume fusion or break-ups 
free parameters chosen to match the MD simulation are P = 1 bar, ji v 
TIP4P water 25J, and po = 0.033 A~ 3 . The coefficient 5 is typically estimated to be between 

291, while VISM was 



0.8 and 1 A for various water models around convex geometries [ll, 
able to predict well the solvation free energies of simple solutes for 8 — 1 A 20j which we 
use in the following. 

We consider ligand positions from d — 11 A to the distance of nearest approach to the 
pocket bottom. The latter is defined as corresponding to a wall-ligand interaction energy 
of 1 ksT and is d ~ —1.8 and -3.8 A for the R5 and R8 system, respectively. We define 
the water occupancy N w of the pocket by the number of oxygens whose LJ centers are 
located at z < 0. Considering the probability distribution P(N W ), we obtain the free energy 
as a function of N w by G(N W ) = — ksT 'In P(N W ) + G'. Without the ligand (effectively 
for d > 9 A), the MD simulation reveals that the R5 pocket is in a stable dry state with 



occupancy N u 



N, 



dry 



0, despite the fact that a few water molecules fit in and consistent 



with experiments on an equally sized protein pocket [6|. The R8 system, however, is found to 
be weakly hydrated. The G(N W ) distribution shown for d = 9 A in Fig. 2 reveals an almost 
barrierless transition between wet and dry states. The metastable wet state comprises 
N w ~ 9 = iV we t water molecules which roughly corresponds to bulk density. 

The approaching ligand considerably changes the G(N W ) distribution in the R8 system. 
As plotted in Fig. 2, for d = 6.5 A the free energy exhibits a minimum at the wet state 
which is slightly stabilized (by ~ 0.4 ksT) over the dry state. The function G(N W ) devel- 
ops, however, concave curvature for N w ~ indicative of the onset of a thermodynamic 



instability. Indeed, upon further approach of the ligand (d = 5.5 A) a local minimum forms 
at the dry state. It becomes a stable, global minimum at the critical distance d c ~ 4.5 A. 
The now metastable wet state completely vanishes for d < A, where we find a free energy 
difference between the wet and dry state of G(iVd r y) — G(N wct ) ~ 5/c^T '. By investigating 
the water density distribution (Fig. 2, right panel), we find that a possible reason for the 
weakly stabilized wet state at d = 6.5 A may be the first methane hydration shell partly 
penetrating the pocket. (This perhaps surprising effect should not be assigned to a lack of 
hydrophobicity or even an hydrophilic nature of the methane but to subtle hydrogen bond 
arrangements within this geometry.) The average occupancy profile (N w (d)) thus exhibits a 
maximum at d — 6.5 A (where (N w ) ~ 6) while it jumps down to ~ at d ~ d c [27]] . 

In the VISM where thermal interface fluctuations are not yet considered, we start the 
numerical relaxation of the functional (1) from (i), one closed solvent boundary which is 
arbitrary and loosely envelopes both the pocketed wall and the ligand, or (ii), the (tight) 
van der Waals surface around the wall and the ligand giving rise to two separated surfaces. 
In Fig. 3 we plot examples of the resulting VISM interfaces for both (i) and (ii), obtained for 
the ligand at d = 4.5 A. For (i) the solution relaxes to a single interface that wraps both wall 
and the ligand together, thereby indicating a dry pocket state [Fig. 3 (a)], while for (ii) the 
solution relaxes to two separate surfaces, one of which closely follows the pocket contours 
indicating a wet state [Fig. 3 (b)]. The existence of two distinct results can be clearly 
attributed to the energy barrier between wet and dry states observed in the simulation (cf. 
Fig. 2). 

By systematically investigating different initial configurations and ligand distances we find 
that the solution for R8 relaxes to at most three distinct interfaces: 1. a single enveloping 
surface around the dry pocket and ligand (Is), 2. two separated surfaces with a dry pocket 
(2s-dry), and 3. two separated surfaces with a wet pocket (2s- wet). Selected examples 
for the interface at ligand distances d = —2,2,4.5 and 9 A are shown in Fig. 3, where we 
plot the bisected VISM interface for a clearer view. For the initial configuration (i) and 
for d < 7 A, the results converge to the Is state while for larger separations a breakup 
into two interfaces (2s-dry) is observed [Fig. 3 (c)]. The stable 2s-dry state exists also 
for 5 < d < 7 A, where it is reached from an initial configuration intermediate between 
(i) and (ii). For the initial configuration (ii) and for d > 0, the results converge to two 
separated surfaces with a wet pocket (2s-wet) while for smaller separations there is only 
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one enveloping surface (Is), see Fig. 3 (d). For R5 we just find two distinct solutions, Is 
and 2s-dry, indicating a very stable dry pocket in agreement with the MD simulations and 
experiments {fj|. These results demonstrate that VISM captures the dewetting transition, 
and the final interface is relaxed into (meta) stable states representing (local) free energy 
minima. This is in physical agreement with the bimodal behavior observed in the MD 
simulation and is further quantified in the following. 

The minimum VISM free energy (1) vs. d is is plotted for R8 in Fig. 4: for d < all 
possible VISM solutions converge to Is, featuring a dry pocket. For < d < d c ~ 4.5 A, 
the 'dry branch' Is is favored over the second appearing branch corresponding to the 2s-wet 
interface (by ~ 8 k^T at d = 0) in excellent agreement with the two-state behavior in the 
MD simulation. For d c < d < 7 A the 2s-dry state is favored over 2s-wet and Is which 
is now highly metastable. For d > 7 A the Is state disappears and 2s-dry is favored by 
roughly 2k bT over 2s- wet. The fact that a dry pocket is favored in VISM for large d is in 
contrast to the MD simulation which supplied a very weakly hydrated pocket for d > 6.5 A. 
Changing the curvature parameter 5 shows that this failure can be attributed to a too high 
energy penalty for concave interface curvature (a too large 6 for H < 0) which favors pocket 
dewetting. It thus appears that the simple curvature correction applied breaks down and is 
not symmetric with respect to concavity and convexity on these small scales. The symmetry 
may be broken by higher order correction terms in the the curvature expansion of the surface 
tension, if feasible [28]. 

If thermal fluctuations were included in VISM, the various energy branches would be 
sampled in a Boltzmann-weighted fashion to yield the solvent-mediated potential of mean 
force (pmf) between the ligand and the pocket. At present, allowing the existence of mul- 
tiple local minima for a given d in the G[v] functional that correspond to the ensemble 
{v } m of most probable solvent configurations, we obtain the ensemble average (EA) as 
G = —ksT In Y2{ v } m e~ G ™ fesT + G", with G" being an arbitrary constant. The resulting 
pmfs G(d) for R8 and R5 are shown in Fig. 4 together with the MD simulation results. The 
curves are overall in good, almost quantitative agreement. A detailed analysis of the indi- 
vidual energy contributions to (1) reveals that the inclusion of the dispersion and curvature 
correction terms in VISM is crucial to capture the onset of the attraction at d c = 4.5 A for 



R8, while SA type of calculations yield a too low d c ~ 



131 ] . Furthermore, the ~ 1 k^T en- 



ergy barrier at d ~ 6 A for the R5 system is nicely captured by VISM; it can be attributed 
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to the unfavorable curvature correction term arising from the development of a concave 
solvent boundary penetrating the pocket, as well as the wall- water dispersion term, whose 
repulsive contribution stems from displacement of water close to the small R5 pocket. An 
EA performed to estimate the average occupancy (N w (d)) for R8 yields qualitative agree- 
ment with the MD, i.e., a step at d ~ d c from zero to nonzero occupancy and a maximum 
at d ~ 6.0 27J]. The step height in VISM (A^ ~ 2), however, is much smaller as in the 



simulation (AN W ~ 6) which is probably due to the continuum approximations of VISM 
and the neglect of fluctuations in performing the EA. 

In summary, the geometry-based VISM is the first implicit solvent model that captures 
the multi-state hydration observed in simulations and experiments [3] and highlights the 



significance of interfacial fluctuations 



291 ] in hydrophobic confinement where the free energy 



can be polymodal. Pocket dewettingmay be regarded as the rate-limiting step for protein- 
ligand binding as found in folding [11]. The existence and height of activation barriers and 
the range of attraction, however, can strongly depend on pocket size and geometry in terms 
of local interface curvature. Our capillary approach, where dispersion and curvature effects 
play explicit roles, may represent a valuable step towards proper interpretation and modeling 
of experimental binding rates. l|. 
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FIG. 1: (a) Sketch of the generic model. The pocket has a radius R. The methane (Me) lig- 
and is fixed at a distance d from the wall surface, (b) MD simulation snapshot illustrating the 
wall/ligand/water system. 
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FIG. 2: MD simulation results for the free energy G(N W ) vs. pocket water occupancy N w in pocket 
R8 for ligand distances d=0.0, 2.5, 4.0, 5.5, 6.5, and 9 A. The curves are shifted vertically and 
we use two scales (1 and 2 /cbT) for a better illustration. The right panel exemplifies the water 
density (p) distribution around pocket and ligand for selected d. 
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FIG. 3: VISM solution of the aqueous interface for the R8 system, a) and b), full three-dimensional 
interface for the ligand at d = 4.5 for one surface (i) and two separated surfaces (ii) as initial 
boundary inputs, respectively, c) and d), the bisected interface for initially one surface (i) and two 
surfaces (ii), respectively, for distances d = -2, 2, 4.5, and 9 A (black, magenta, blue, and red). 
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FIG. 4: VISM free energies for the Is (squares), 2s-wet (circles), and 2s-dry (triangles) branches for 
R8 (top) and R5 (bottom) and the solvent-mediated pmf between the pocket and ligand from MD 
simulations (solid lines) and the ensemble average (EA) over the VISM branches (filled circles). 
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The wall-water interaction 



In order to construct hydrophobic walls we considered 
a paraffine-like material of 0.8 g/cm 3 density composed 
of CH 2 units. Assuming a hexagonally closed-packed ar- 
rangement, the given density requires a lattice constant 
of 3.5 A which is too coarse to produce a relatively 
smooth hemispherical pocket. Thus, we reduced the lat- 
tice constant to 1.25 A while at the same time adjust- 
ing the Lennard- Jones potential parameters of the wall 
pseudo-atoms to reproduce the original parafhne wall - 
water interaction energy (see inset to Fig. [I]) that was ob- 
tained with the united atom OPLS parameters for CH2 
units [1]. The wall- water interaction energy was calcu- 
lated by simply averaging the interaction energy over 
planes at constant z. 




z/A 

FIG. 1: Water oxygen density vs. the distance z from the 
flat wall surface; g(z) = 1.0 corresponds to a water number 
density of 0.0334 A -3 . Inset: wall- water interaction energy 
for the original 3.5 A grid lattice (dashed line), and the 1.25 A 
grid lattice with adjusted LJ potential parameters (squares). 



Average water occupancy of the pocket 

Based on the VISM results we estimated an average 
pocket occupancy (N w ) for different ligand positions by 
the ensemble average 



(N w ) = N[v]e- G[v]/kaT / 



-G[v]/k B T 



{v} m 



(1) 



where N[v] = for dry- type solutions and N[v] = 
N we t = 9 for wet-type solutions. A comparison with 
MD results (Fig. [2]) shows the correct qualitative behav- 
ior with a transient increase in pocket wetting at d ~ 6 A, 
followed by ligand induced dewetting. The quantitative 
discrepancy is probably due to a) the approximations in 
the VISM functional leading to over-stabilization of dry 
state relative to wet state, and b) including only local 
G[v] minima in the ensemble average while omitting in- 
termediate states of not much higher free energy. We 
expect the qualitative agreement to improve upon inclu- 
sion of interface thermal fluctuations into VISM. 
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FIG. 2: Average pocket occupancy (N w ) from MD simulation 
and VISM ensemble average. 



The height of the first peak in the wall-water density 
profile from our MD simulations (Fig. [T]) is within the 
range (1.3 to 1.6) observed in all atom MD simulations 
of hydrocarbon- water interfaces [ESS], suggesting that 

the walls indeed closely resemble a paraffine-like material. * Electronic address: jdzubiel@ph.tum.de 
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